Nathan Trouvain
Inria - Mnemosyne
PFIA - 28 juin 2021
list[sequence], ou une séquence représente une série temporelle.Série temporelle de Mackey-Glass
L'équation de Mackey-Glass décrit l'évolution de différent phénomènes biologiques, comme la quantité relative de globules rouges dans le temps.
Cette équation est une équation différentielle retardée :
$$ \frac{dP(t)}{dt} = \frac{a P(t - \tau)}{1 + P(t - \tau)^n} - bP(t) $$où $a = 0.2$, $b = 0.1$, $n = 10$. Le retard $\tau$ est de $17$ pas de temps. $\tau$ a une forte influence sur le caractère chaotique de la série temporelle modelée. $17$ donne des résultats plutôt chaotiques.
from reservoirpy.datasets import mackey_glass
timesteps = 2510
tau = 17
X = mackey_glass(timesteps, tau=tau)
# rescale between -1 and 1
X = 2 * (X - X.min()) / (X.max() - X.min()) - 1
plot_mackey_glass(X, 500, tau)
Prédire $P(t + 10)$ connaissant $P(t)$.
from reservoirpy.datasets import to_forecasting
x, y = to_forecasting(X, forecast=10)
X_train1, y_train1 = x[:2000], y[:2000]
X_test1, y_test1 = x[2000:], y[2000:]
plot_train_test(X_train1, y_train1, X_test1, y_test1)
units = 100
leak_rate = 0.3
spectral_radius = 1.25
input_scaling = 1.0
connectivity = 0.1
input_connectivity = 0.2
regularization = 1e-8
seed = 1234
Win = mat_gen.generate_input_weights(units, 1, input_scaling=input_scaling,
proba=input_connectivity, input_bias=True,
seed=seed)
W = mat_gen.generate_internal_weights(units, sr=spectral_radius,
proba=connectivity, seed=seed)
reservoir = ESN(leak_rate, W, Win, ridge=regularization)
L'apprentissage est offline ("hors-ligne") : il n'a lieu qu'une seule fois, sur l'ensemble des données d'entraînement.
reservoir.train([X_train1], [y_train1], verbose=True)
Training on 1 inputs (2000 steps) -- wash: 0 steps
Train: 100%|██████████| 2000/2000 [00:00<00:00, 2537.93it/s]
reservoir.Wout is not None
True
y_pred1, states1 = reservoir.run([X_test1], verbose=True)
y_pred1 = y_pred1[0]
Running on 1 inputs (500 steps)
Run: 100%|██████████| 500/500 [00:00<00:00, 1439.26it/s]
plot_results(y_pred1, y_test1)
Coefficient de détermination $R^2$ et erreur quadratique normalisée :
r2_score(y_test1, y_pred1), nrmse(y_test1, y_pred1)
(0.444170605212862, 0.23219740216886683)
states = reservoir.train([X_train1], [y_train1], return_states=True)
last_state = states[0][-1]
y_pred1, _ = reservoir.run([X_test1], init_state=last_state)
y_pred1 = y_pred1[0]
plot_results(y_pred1, y_test1)
Coefficient de détermination $R^2$ et erreur quadratique normalisée :
r2_score(y_test1, y_pred1), nrmse(y_test1, y_pred1)
(0.9999562183846586, 0.003153349937792222)
Passons d'un horizon de prédiction de 10 pas de temps à un horizon de 100 pas de temps
x, y = to_forecasting(X, forecast=100)
X_train2, y_train2 = x[:2000], y[:2000]
X_test2, y_test2 = x[2000:], y[2000:]
plot_train_test(X_train2, y_train2, X_test2, y_test2)
plot_results(y_pred2, y_test2, sample=400)
Determination coefficient $R^2$ and NRMSE:
r2_score(y_test2, y_pred2), nrmse(y_test2, y_pred2)
(0.9686082619494742, 0.06591305950469496)
units = 500
leak_rate = 0.3 # - leaking rate
spectral_radius = 0.99 # - rayon spectral
input_scaling = 1.0 # - facteur de mise à l'échelle des entrées (input scaling)
connectivity = 0.1 # - densité des connexions du reservoir vers lui même
input_connectivity = 0.2 # et des entrées vers le reservoir
regularization = 1e-4 # - coefficient de régularisation (L2)
seed = 1234 # reproductibilité
reservoir = reset_esn()
x, y = to_forecasting(X, forecast=1)
X_train3, y_train3 = x[:2000], y[:2000]
X_test3, y_test3 = x[2000:], y[2000:]
reservoir.train([X_train3], [y_train3])
start = 0; seed_timesteps = 100; nb_generations = 400
warming_inputs = X_test3[start:start+seed_timesteps]
X_t = X_test3[start+seed_timesteps: start+nb_generations+seed_timesteps]
X_gen, states, warming_out, warming_states = reservoir.generate(nb_generations, warming_inputs)
plot_generation(X_gen, X_t, nb_generations, warming_out=warming_out,
warming_inputs=warming_inputs, seed_timesteps=seed_timesteps)
states = reservoir.train([X_train3], [y_train3], return_states=True)
X_gen, states, _, _ = reservoir.generate(nb_generations, init_state=states[0][-1])
plot_generation(X_gen, X_test3[:nb_generations], nb_generations)
Apprentissage se déroulant de manière incrémentale.
Utilisation de l'algorithme FORCE (Sussillo and Abott, 2009)
from reservoirpy import ESNOnline
Win = mat_gen.generate_input_weights(units, 1, input_scaling=input_scaling,
proba=input_connectivity, input_bias=True,
seed=seed)
W = mat_gen.generate_internal_weights(units, sr=spectral_radius,
proba=connectivity, seed=seed)
reservoir_online = ESNOnline(leak_rate, W, Win, dim_out=1)
reservoir_online.Wout is not None # Wout est créée dès le début !
True
Le dernier état calculé est cette fois-ci stocké dans l'objet. (Parallélisation impossible)
reservoir_online.state is not None
True
outputs_pre = np.zeros(X_train1.shape)
for t, (x, y) in enumerate(zip(X_train1, y_train1)): # pour chaque pas de temps de la série :
# sortie avant entraînement, mise à jour de reservoir.state
output_pre, state = reservoir_online.compute_output(x.reshape(1, -1))
# entraînement à partir de reservoir.state
reservoir_online.train_from_current_state(targeted_output=y.reshape(1, -1))
outputs_pre[t, :] = output_pre
plot_results(outputs_pre, y_train1, sample=100)
plot_results(outputs_pre, y_train1, sample=500)
reservoir_online.reset_reservoir()
states = reservoir_online.train([X_train1], [y_train1])
pred_online, _ = reservoir_online.run([X_test1]) # Wout est maintenant figée
pred_online = pred_online[0]
plot_results(pred_online, y_test1, sample=500)
Determination coefficient $R^2$ and NRMSE:
r2_score(y_test1, pred_online), nrmse(y_test1, pred_online)
(0.9999602106090204, 0.006511751721725012)
Les données peuvent être téléchargées sur Zenodo : https://zenodo.org/record/4736597
im = plt.imread("./static/canary.png")
plt.figure(figsize=(5, 5)); plt.imshow(im); plt.axis('off'); plt.show()
display(audio)
Plusieurs motifs temporels répétitifs differents à décoder : les phrases.
im = plt.imread("./static/canary_outputs.png")
plt.figure(figsize=(15, 15)); plt.imshow(im); plt.axis('off'); plt.show()
X_train, y_train = X[:-10], Y[:-10]
X_test, y_test = X[-10:], Y[-10:]
reservoir.train(X_train, y_train, verbose=True, workers=-1)
Training on 90 inputs (266054 steps) -- wash: 0 steps
Train: 100%|██████████| 266054/266054 [00:23<00:00, 11121.76it/s]
outputs, _ = reservoir.run(X_test, verbose=True)
Running on 10 inputs (32552 steps)
Run: 100%|██████████| 32552/32552 [00:02<00:00, 11873.05it/s]
scores # pour chaque chant testé
[0.9370629370629371, 0.9093011305241521, 0.9499136442141624, 0.9412306983175847, 0.962171052631579, 0.9540372670807453, 0.8858530661809351, 0.9371994342291372, 0.9174974217944311, 0.8932743921692453]
print("Précision moyenne :", f"{np.mean(scores):.4f}", "±", f"{np.std(scores):.5f}")
Précision moyenne : 0.9288 ± 0.02471
units = 100 # - nombre de neurones dans le reservoir
leak_rate = 0.3 # - leaking rate
spectral_radius = 1.25 # - rayon spectral
input_scaling = 1.0 # - facteur de mise à l'échelle des entrées
connectivity = 0.1 # - densité des connexions du reservoir vers lui même
input_connectivity = 0.2 # et des entrées vers le reservoir
regularization = 1e-8 # - coefficient de régularisation (L2)
seed = 1234 # reproductibilité
Le rayon spectral est la valeur propre maximale de la matrice des poids du réservoir ($W$).
reservoir = reset_esn()
states = []
radii = [0.1, 1.25, 10.0]
for sr in radii:
Win = mat_gen.generate_input_weights(units, 1, input_scaling=0.1, input_bias=True)
W = mat_gen.generate_internal_weights(units, sr=sr,
proba=connectivity, seed=seed)
reservoir.W = W
reservoir.Win = Win
s = reservoir.compute_all_states([X_test1[:500]])
states.append(s[0])
units_nb = 20
plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
plt.subplot(len(radii)*100+10+i+1)
plt.plot(s[:, :units_nb], alpha=0.6)
plt.ylabel(f"$sr={radii[i]}$")
plt.xlabel(f"Activations ({units_nb} neurons)")
plt.show()
$-$ rayon spectral $\rightarrow$ dynamiques stables
$+$ rayon spectral $\rightarrow$ dynamiques chaotiques
Rayon spectral et Echo State Property : rayon spectral $\rightarrow$ 1 (assure que les états internes ne sont pas affectés par l'initialisation).
Il s'agit d'un coefficient appliqué à $W_{in}$, venant changer l'échelle des données en entrée.
reservoir = reset_esn()
states = []
scalings = [0.01, 0.1, 1.]
for iss in scalings:
Win = mat_gen.generate_input_weights(units, 1, input_scaling=iss,
proba=input_connectivity, input_bias=True,
seed=seed)
reservoir.Win = Win
s = reservoir.compute_all_states([X_test1[:500]])
states.append(s[0])
units_nb = 20
plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
plt.subplot(len(scalings)*100+10+i+1)
plt.plot(s[:, :units_nb], alpha=0.6)
plt.ylabel(f"$iss={scalings[i]}$")
plt.xlabel(f"Activations ({units_nb} neurons)")
plt.show()
Correlation moyenne des activités des neurones du reservoir avec les entrées :
for i, s in enumerate(states):
corr = correlation(states[i], X_test1[:500])
print(f"ISS : {scalings[i]}, correlation moyenne : {corr}")
ISS : 0.01, correlation moyenne : -0.032300266438732315 ISS : 0.1, correlation moyenne : 0.5455501730434691 ISS : 1.0, correlation moyenne : 2.441693673875052
L'input scaling peut aussi être utilisé pour ajuster l'influence de chaque donnée en entrée.
avec $\alpha \in [0, 1]$ et:
$$ f(u, x) = \tanh(W_{in} \cdotp u + W \cdotp x) $$reservoir = reset_esn()
Win = mat_gen.generate_input_weights(units, 1, input_scaling=0.1, input_bias=True)
reservoir.Win = Win
states = []
rates = [0.02, 0.2, 0.9]
for lr in rates:
reservoir.lr = lr
s = reservoir.compute_all_states([X_test1[:500]])
states.append(s[0])
units_nb = 20
plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
plt.subplot(len(rates)*100+10+i+1)
plt.plot(s[:, :units_nb] + 2*i)
plt.ylabel(f"$lr={rates[i]}$")
plt.xlabel(f"States ({units_nb} neurons)")
plt.show()
Le leaking rate contrôle la "mémoire" de l'ESN. Il peut être vu comme l'inverse de sa constante de temps